%% SIMULATION PERSONALIZED DEFAULT SETTING
clc;
clear all;

%% SECTIONS
% 1. set default= rho + a
% 2. set default= rho * b

%% Section 1: set default= rho + a
data=csvread('data_section4.csv',1,0);

% use lambdas obtained by bootstrapping, N=1000
load('lambdas_bootstrapping_1000_global.mat')

% use subset of AD treatment for optimization >0
default=data(:,1);
donation=data(:,2);
logInd = default== 0;
rho_AD=donation(logInd);
rho_AD=rho_AD(rho_AD>0);
rho_AD=round(rho_AD);
rho_AD(rho_AD==0)=1;

%% set seed
rng(432123486);
%% take random sample from AD treatments
N=100000;
rho = randsample(rho_AD,N, true);
%% Estimating a*
% for each estimated set of lambdas obtained by bootstrapping, find a* that
% maximizes total donations.
c=1000;
as=zeros(c,1);
totals=zeros(c,1);

for k=1:c
lambda1=lambdas(k,1);
lambda2=lambdas(k,2);
lambda3=lambdas(k,3);
indiv_don=@(a)negtotaldonationsadd(rho, lambda1, lambda2,lambda3, a);
opts=optimoptions(@fmincon,'Algorithm','interior-point', 'FunctionTolerance',1e-12, 'StepTolerance',1e-8,'AlwaysHonorConstraints','bounds', 'MaxFunEvals',2000);
problem = createOptimProblem('fmincon','x0', 2 , 'objective', indiv_don, ...
                             'lb', 0, 'ub', Inf,'options', opts);
gs = GlobalSearch;
[as(k,1),totals(k,1)] = run(gs,problem);
end
a_star=mean(as);
total_add=mean(totals);

display('to optimize the total donation through adding an amount, a, to rho; set a:');
disp(a_star');
     
display('total simulated donations at default=rho+a*')
disp(-total_add);

display('total donation at default=rho')
disp(sum(rho));


add= [a_star, -total_add, sum(rho)]';
T=table(add);
writetable(T,'_results/results_pers_defaults_optim_bootstrap_add.xls')

%% Figure Donations under Personalized Defaults Additive
clc;
clear all;

data=csvread('data_section4.csv',1,0);
% use lambdas obtained from bootstrapping to construct the confidence
% interval
load('lambdas_bootstrapping_1000_global.mat')

% use all observations of AD
default=data(:,1);
donation=data(:,2);
logInd = default== 0;
rho_AD=donation(logInd);
rho_AD=round(rho_AD);

%% set seed
rng(4321656);
%% take relfreq as underlying distribution: simulate rhos accordingly
% use large sample for simulation to approximate well the alpha and delta distributions
%The sample size is set similar to the sample size used for the simulations of 
% donations with 10, 20, and 50 defaults where N=500.000, 
%but since now we use all observations including the zero donations 
% we rescale N accordingly: N= 500.000*(#total observations / #non-zero observations).


p=length(rho_AD(rho_AD>0));
c= 500000/p;
N=c*length(rho_AD);
N=round(N);
rho = randsample(rho_AD,N, true);
mean(rho);
rhos=rho(rho>0);
addz=length(rho)-length(rhos);
%%
aa=80
means=zeros(aa,1000);
total=zeros(aa,1000);
for k=1:1000;
 lambda1=lambdas(k,1);
 lambda2=lambdas(k,2);
 lambda3=lambdas(k,3);

n=length(rhos);
z=zeros(n,1);
z=[rand(1,n)>lambda1];
nd=sum(z);
z=(z==1);

delta=zeros(n,1);
mu=1/lambda2;
delta(z) = exprnd(mu, [nd, 1]); %delta exponentially distributed with mean 1/lambda2

alpha=zeros(n,1);
mu3=1/lambda3;
alpha(z) = exprnd(mu3, [nd, 1]);
for i=1:aa;
  [means(i,k) total(i,k)]=default_setting_add_wz(rhos, alpha, delta,addz, z, n, i);
end
end

meandon=zeros(aa,1);
SEM=zeros(aa,1);
CI=zeros(aa,2);
ts=tinv([0.025  0.975],99);
for i=1:aa;
meandon(i,1)=mean(means(i,:));
SEM(i,1) = std(means(i,:))/sqrt(length(1000)); 
CI(i,1:2) = mean(means(i,:)) + ts*SEM(i,1);
end 
save('personalizing_default_add_wz_results.mat', 'means', 'meandon', 'CI','lambdas');

%% safe figure
load('personalizing_default_add_wz_results.mat');
figure
X=[1:80];
plot(X,meandon,'k','LineWidth',1)
hold on
ciplot_trans(CI(:,1),CI(:,2),X, 'k--', 0);
set(gca,'XLim',[1 num]);
set(gca,'YLim',[1.4 2]);
set(gca,'fontsize',18)
ya=ylabel('Mean donation in Euros','fontsize',22)
xa=xlabel('\it a','Interpreter','tex','fontsize',26)
hold off
h=gcf;
set(h,'PaperOrientation','landscape');
set(h,'PaperUnits','normalized');
set(h,'PaperPosition', [0 0 1 1]);
print(gcf, '-dpdf', 'individualized_default_setting_add_bootstrap_wz.pdf');

%% SECTION 2 %% set default=rho*b
% for each estimated set of lambdas from the bootstrapping, find b* that
% maximizes total donations.     
clc;
clear all;

data=csvread('data_section4.csv',1,0);
load('lambdas_bootstrapping_1000_global.mat')
% use subset of AD treatment for optimization >0
default=data(:,1);
donation=data(:,2);
logInd = default== 0;
rho_AD=donation(logInd);
rho_AD=rho_AD(rho_AD>0);
rho_AD=round(rho_AD);
rho_AD(rho_AD==0)=1;
     
%% set seed
rng(432856487);
%% take random sample from AD treatments
N=100000;
rho = randsample(rho_AD,N, true);

c=1000;
ps=zeros(c,1);
totals=zeros(c,1);
     
for k=1:c;
lambda1=lambdas(k,1);
lambda2=lambdas(k,2);
lambda3=lambdas(k,3);
indiv_don=@(b)negtotaldonationsmult(rho, lambda1, lambda2,lambda3, b);
opts=optimoptions(@fmincon,'Algorithm','interior-point', 'FunctionTolerance',1e-12, 'StepTolerance',1e-8,'AlwaysHonorConstraints','bounds', 'MaxFunEvals',2000);
problem = createOptimProblem('fmincon','x0', 2 , 'objective', indiv_don, ...
                                  'lb', 0, 'ub', Inf,'options', opts);
gs = GlobalSearch;
[bs(k,1),totals(k,1)] = run(gs,problem);
end
     
b_star=mean(bs);
total_mult=mean(totals);
display('to optimize the total donation through multiplying rho with b; set b:');
disp(b_star');
          
display('total simulated donations are ')
disp(-total_mult);
          
display('total donation in the AD treatment are')
disp(sum(rho));
    
mult= [b_star, -total_mult, sum(rho)]';
T=table(mult);
writetable(T,'_results/results_pers_defaults_optim_bootstrap_mult.xls')
     
%% Figure Donations under Personalized Defaults Multiplicative
clc;
clear all;
     
     
data=csvread('data_section4.csv',1,0);
     
load('lambdas_bootstrapping_1000_global.mat')
     
% use all observations AD
default=data(:,1);
donation=data(:,2);
logInd = default== 0;
rho_AD=donation(logInd);
rho_AD=round(rho_AD);
     

%% take relfreq as underlying distribution: simulate rhos accordingly
% use large sample for simulation to approximate well the alhpa and delta distributions
% similar sample size as the simulations with 10, 20, and 50 default where we used N=500000, 
% but since now we use all observations including the zero donations:
% N= 500.000*(#total observations / # non-zero observations)

p=length(rho_AD(rho_AD>0));
c= 500000/p;
N=c*length(rho_AD);
N=round(N);
rho = randsample(rho_AD,N, true);
mean(rho);
rhos=rho(rho>0);
addz=length(rho)-length(rhos);
a1=(1:0.01:4);
bb=length(a1)
means=zeros(bb,1000);
total=zeros(bb,1000);

for k=1:1000;
 lambda1=lambdas(k,1);
 lambda2=lambdas(k,2);
 lambda3=lambdas(k,3);

n=length(rhos);
z=zeros(n,1);
z=[rand(1,n)>lambda1];
nd=sum(z);
z=(z==1);

delta=zeros(n,1);
mu=1/lambda2;
delta(z) = exprnd(mu, [nd, 1]); %delta exponentially distributed with mean 1/lambda2

alpha=zeros(n,1);
mu3=1/lambda3;
alpha(z) = exprnd(mu3, [nd, 1]);
for i= 1:bb;
  [means(i,k) total(i,k)]=default_setting_mult_wz(rhos, alpha, delta, addz, z, n, a1(i));
end
end

meandon=zeros(bb,1);
SEM=zeros(bb,1);
CI=zeros(bb,2);
ts=tinv([0.025  0.975],99);
for i=1:bb;
meandon(i,1)=mean(means(i,:));
SEM(i,1) = std(means(i,:))/sqrt(length(1000)); 
CI(i,1:2) = mean(means(i,:)) + ts*SEM(i,1);
end 

save('personalizing_default_mult__wz_results.mat', 'means', 'meandon', 'CI','lambdas');

%% safe figure

load('personalizing_default_mult__wz_results.mat')
figure
plot(a1,meandon,'k','LineWidth',1)
hold on
ciplot_trans(CI(:,1),CI(:,2),a1, 'k--', 0);
set(gca,'YLim',[1.4 2]);
set(gca,'fontsize',18)
ya=ylabel('Mean donation in Euros','fontsize',22)
xa=xlabel('\it b','Interpreter','tex','fontsize',26)
hold off
h=gcf;
set(h,'PaperOrientation','landscape');
set(h,'PaperUnits','normalized');
set(h,'PaperPosition', [0 0 1 1]);
print(gcf, '-dpdf', 'individualized_default_setting_mult_bootstrap_wz.pdf');



